Study of first-order interface localization-delocalization transition 
in thin Ising-films using Wang-Landau sampling 
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Using extensive Monte Carlo simulations, we study the interface localization-delocalization tran- 
sition of a thin Ising film with antisymmetric competing walls for a set of parameters where the 
transition is strongly first-order. This is achieved by estimating the density of states (DOS) of the 
model by means of Wang-Landau sampling (WLS) in the space of energy, using both, single-spin- 
flip as well as N-fold way updates. From the DOS we calculate canonical averages related to the 
configurational energy, like the internal energy, the specific heat, as well as the free energy and the 
entropy. By sampling microcanonical averages during simulations we also compute thermodynamic 
quantities related to magnetization like the reduced fourth order cumulant of the order parameter. 
We estimate the triple temperatures of infinitely large systems for three different film thicknesses 
via finite size scaling of the positions of the maxima of the specific heat, the minima of the cumulant 
and the equal weight criterion for the energy probability distribution. The wetting temperature 
of the semi-infinite system is computed with help of the Young equation. In the limit of large 
film thicknesses the triple temperatures are seen to converge towards the wetting temperature of 
the corresponding semi-infinite Ising model in accordance with standard capillary wave theory. We 
discuss the slowing down of WLS in energy space as observed for the larger film thicknesses and 
lateral linear dimensions. In case of WLS in the space of total magnetization we find evidence that 
the slowing down is reduced and can be attributed to persisting free energy barriers due to shape 
transitions. 



INTRODUCTION 



The restriction of the geometry of a condensed-matter 
system has fundamental impact on a phase transition. In 
a finite system, sharp phase transitions can no longer oc- 
cur, since the free energy is then an analytic function of 
its independent variables and the transition is rounded off 
and shifted. A particular realization of a confined geom- 
etry in d = 3 dimensions, playing a pivotal role due to its 
fundamental importance in material science and technol- 
ogy, are thin films, infinitely extended in two directions 
but of finite thickness D, where the transition is now 
not only shifted away from its bulk value, corresponding 
to D — > do, but also changes its character from three- 
dimensional to two-dimensional. As an example we may 
consider here a fluid near a gas-liquid coexistence in the 
bulk, or similarly, an (A,B) binary mixture or alloy near 
phase coexistence, confined between two parallel walls. 
Of particular interest is the case, where the two walls of 
the system prefer different phases, i.e., one wall favors 
high-density liquid (or A-particles) while the other one 
prefers low-density gas (or B-particles) , which is com- 
monly termed "competing walls" . A generic model for 
such systems actually is the nearest neighbor Ising model 
in a thin film geometry where one now has two surfaces 
a distance D apart, on which magnetic surface fields 
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Hi = —Hd of opposite sign but equal magnitude act in 
order to mimic the competing walls (cf. Fig.^l. In addi- 
tion one allows for a different interaction J s > between 
nearest neighbors located in the surfaces, while nearest 
neighbors in the bulk interact with a coupling J > 0. 
The meaning of the magnetic surface fields becomes ap- 
parent, when reinterpreting the Ising Hamiltonian as a 
lattice gas for a fluid or a mixture, where Ising spins 
Si = — 1 or Si = +1 now correspond to lattice sites i be- 
ing empty or occupied, or being taken by an A-particle or 
a i?-particle, respectively. Then, surface magnetic fields 
translate into chemical potentials, i.e., binding energies 
to the walls. 

Remarkably, the transition that one encounters in the 
Ising film differs from the transition in the bulk system 
at T cb 0, H H S B 0, H : For all finite thicknesses D 
of the film, the transition at T c b is completely rounded 
off and no singular behavior shows up, despite the fact 
that the system is infinite in the other directions. In- 
stead, one observes a transition at a lower temperature 
Tq(D) < T c b, at which the system changes from a state 
with a delocalized interface running parallel to the walls 
in the center of the film (T (D) < T < T ch ), to a twofold 
degenerate state (T < T (D)), where the interface is now 
localized near one of the two walls (cf . Fig. ^ . 
Most interestingly, for D — > oo, the transition temper- 
ature To(D) of the interface localization-delocalization 
does not converge towards the bulk critical temperature 
T^b, but towards the wetting temperature T w (Hi) at 
which a macroscopically thick liquid layer (spins point- 
ing upwards) wets the surface in the corresponding semi- 
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FIG. 1: (a) Thin film geometry with two free surfaces at n — 1 
and n = D (shaded gray) on which magnetic surface fields Hi 
and Hd act. Here, the surface at n = 1 favors spin up (+), 
while the surface at n = D favors spin down (-). Parallel to 
the LxL surfaces, periodic boundary conditions are imposed, 
(b) Delocalized Interface, (c) Interface located at either of the 
two surfaces. 



infinite system. Thus, the nature of the transition at 
finite D is seen to depend on the nature of the wetting 
transition in the underlying semi-infinite system. Upon 
enhancing the interaction J s of spins in the surfaces with 
respect to the bulk interaction J one can tune the wetting 
transition and thus the interface transition for finite film 
thicknesses D to be of first order Q, i.e., T (D) = T tI (D) 
is now a triple point where the three phases shown in 
Fig. [IJb)-(c) coexist. By reducing the film thickness one 
may then pass through a tricritical point where the or- 
der of the transition changes from first to second order 

Our paper is arranged as follows: First, we briefly in- 
troduce the thin film Hamiltonian and give a description 
of the employed Wang-Landau sampling (WLS) which 
aims at sampling the density of states (DOS) directly. 
The slowing down of WLS for our model, encountered 
especially for large system sizes is discussed. With re- 
gard to these difficulties we then propose to split the 
DOS in a branch contributing to the ordered phase and 
one contributing to the disordered phase, which we nor- 
malize separately. We then present the thermodynamic 
quantities calculated from the DOS and compute the in- 
finite lattice triple temperatures from the various finite 
size data. Finally, the wetting temperature of the semi- 
infinite system is determined via the Young equation and 
the convergence of the triple temperatures towards the 
wetting temperature for increasing film thickness is ex- 
amined. We close with a brief discussion of our results. 



II. MODEL AND SIMULATION METHOD 

We consider the Ising Hamiltonian on a cubic lattice 
in a L x L x D geometry (cf. Fig.^a)), where N = L 2 D 
is the total number of spins Sf. 

(i,j)h (i,j)s i 



ifEsurfacc 1 



zfEsurfacc D 



Here, the sum (i,j)b runs over all pairs of nearest neigh- 
bors where at least one spin is not located in one of the 
surfaces and the sum (i,j) s runs over all pairs of near- 
est neighbors with both spins located in one of the two 
surfaces. In this paper we study three different film thick- 
nesses D = 6, 8, 12, and linear lateral dimensions ranging 
from L — 16 to L = 128 (for the two largest choices of 
D the minimal L is L = 32 and L = 48, respectively). 
We restrict ourselves here to antisymmetric surface fields 
Hi = —Hd and bulk field, H — 0. By virtue of the sym- 
metry there is an exact degeneracy of the phases where 
the interface is bound to either of the surfaces, and the 
triple point and the phase coexistence below Tq(D) oc- 
curs at H = 0. We do not study pre- wetting like phase 
coexistence for T > Tq(D) and H ^ 0. Specifically we 
choose Hi/ J — 0.25 and J s /J = 1.5. For these parame- 
ters the interface localization-delocalization transition is 
clearly first-order for all thicknesses D. Already for a 
smaller surface-to-bulk coupling ratio J s /J — 1.45, the 
transition turned out to be so strongly first order accord- 
ing to the study of Ferrenberg et. al 0, that lattices 
with D = 8 and L > 32 could not be equilibrated using 
a standard canonical heatbath algorithm. The reason 
for such difficulties can be seen directly from the canoni- 
cal probability distribution, Pl,d(E), of the energy that 
develops two pronounced peaks at the transition point, 
corresponding to the coexisting ordered (— ) and disor- 
dered phases (+) which are separated by a deep min- 
imum P^ L1 ^ ) (E) corresponding to the mixed phase con- 
figurations (cf. FigJ2J. Here, one has additional in- 
terfaces in the system which cost an extra free energy 
AFl d — "fDL, where 7 is of the order of the inter- 
face tension between the two oppositely oriented domains 
of spins. This yields Pjj?$(E) oc exp(-/5A^ L , D ), where 
[3 = denotes the inverse temperature. Hence, any 

simulation technique which aims at sampling a canoni- 
cal energy probability distribution oc g(E) exp(— E/ksT) 
directly will become trapped in the phase in which the 
system was initially prepared and may practically never 
escape from there, even in case of relatively small sys- 
tems. 

In order to give an example for the strong metastabil- 
ity, Fig. |3 shows hysteresis-loops of the internal energy 
per spin (e) = (E)/N which were recorded using a con- 
ventional Metropolis Monte Carlo algorithm for a sys- 
tem of size D — 12 and L — 48. The simulations were 
started in the disordered phase. In case cooling is per- 
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— D=6; L=64, 96, 128 

— D=8; L=64, 96, 128 " 




FIG. 2: Energy probability distributions Pl,d at equal 
weight. The peak positions e~(L, D) and e + (L, D) (indicated 
for D = 6 and L = 128) define the finite volume latent heats 
Ae(L,D) = e + (L,D) — e~{L,D). Arrows pointing on the 
energy axis indicate the interval 7 CO nter, Eq. HUH , in case of 
L = 128 and D = 6. 

formed too fast (open circles in Fig. [3J| one reaches the 
roughening temperature Tr while still being in the dis- 
ordered soft-mode phase, i.e., the interface becomes flat 
in the center of the film and it becomes impossible to 
reach the ordered phase upon further cooling. Using a 
much larger simulational effort (~ 10 7 MCS) one obtains 
a closed loop - although the observed hysteresis is still 
huge - which clearly indicates a phase transition in the 
range 0.244 < J/3 tr < 0.341. Locating the exact transi- 
tion point in this way would however require an enormous 
simulational effort even for the moderate system size at 
hand. An improvement results from thermodynamic in- 
tegration of the low- and high-temperature branches of 
the internal energy, which yields the free energy per site 
/ Hi: 

0f(J3) = ft ef /(/3ref) + / (e^/d/?. (2) 

J Pre! 

For the reference values we have regarded the spins 
as noninteracting at J/3 ro f = 0.00005, i.e., /(/3 re f) = 
— P~] In 2, while on the low temperature side the free en- 
ergy was matched with a series expansion based on the 
first two excited states at J/3 rc f = 1.10005. The crossing 
point of both branches of the free energy then yields the 
transition point, which can be determined with an accu- 
racy of 0.4%. 

The result, that the correct location of the first order 
transition is not in the middle of the hysteresis loop but 
very close to its end at the high temperature side (dashed 
curve in Fig. is very surprising at first sight. It should 
be noted however, that the hysteresis observed in Monte 
Carlo simulations has nothing to do with the "Maxwell 



equal area rule" of mean field theories, but is of kinetic 
origin. The almost free interface in the center of the film 
is very slowly relaxing and feels only a very weak po- 
tential from the walls, and thus is much more metastable 
than the state where the interface is tightly bound to one 
of the walls. 



A. Wang-Landau Sampling 

In order to avoid problems due to metastability and to 
further increase accuracy, we have decided to use Wang- 
Landau sampling (WLS) dElEE! 

in order to com- 
pute thermodynamic quantities of the systems via esti- 
mating the density of states (DOS) of Hamiltonian Q. 
In WLS one accepts trial configurations with probability 
mm[l, g(E)/g(E')}, where g(E) is the DOS and E and 
E' are the energies of the current and the proposed con- 
figuration, respectively. At each spin flip trial the DOS 
is modified g(E) — * g(E) ■ fi by means of a modification 
factor fi, which is reduced according to fi + \ = y/Jl in 
case the recorded energy histogram H{E) is flat within 
some percentage e of the average energy histogram, i.e., 
H{E) > e{H(E')) E > for all E. H(E) is then reset to 
zero and the procedure is repeated until a flat H(E) is 
achieved using a final modification factor /fi na i. In prac- 
tice one samples a logarithm of the DOS, i.e., \og w g(E) 
since g{E) may become very large and modifying the 
DOS then corresponds to adding a small modification 
increment Asi — log 10 /$ . The implementation of the 
single-spin-flip WLS is straightforward and we refer the 
reader to Refs. |l3lll4l ll5| for details. When considering 
systems with a large number of distinct energy levels it is 
useful to partition the entire energy range into adjacent 
subintervals in order to sample the DOS in a parallel 
fashion. For energy intervals that contain states with 
low degeneracy, e.g., the ground state, one can further 
accelerate WLS by combining it with the rejection-free 
N-fold way of Bortz et. al [la. Il7|. Here, the underly- 
ing idea is to partition all spins S v , v £ {1, ■■•,N} into 
M classes c„ G {0, M — 1} according to the change in 
energy AE Cu caused by flipping a spin S v at site v. Mak- 
ing explicit use of H\ = —Hp and H = in the Ising 
Hamiltonian QJ, we can evaluate AE Cv as follows: 

!2(Ju„ + J s v v + Hi)S u if v e surfc.l 
2(Ju v + J s v v -Hx)S v if vesmic.D (3) 
2Ju v S v else, 

where S v is the spin value before it is overturned and u v 
and v v denote sums over the nearest neighbor sites niy) 
of site v. 

!J2 ^ v <£ surface 

J2 S^) if v e surface, ^ 
\x{y) ^surface 
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FIG. 3: (a) Energy hysteresis curves. Cooling (heating) was performed at a rate of JA/3/MCS = 4.3403 ■ 10 6 (open circles, 
note that not all data points are plotted) and in steps of JA/3 = 0.0005, using 100 MCS for equilibration at each j3 and another 
10 4 MCS for measuring the energy (solid line) . The equilibrium curve obtained from WLS in the space of energy is also shown. 
The roughening temperature J/k&Tn — 0.40758(1) flOl) is indicated by an arrow, (b) Low and high temperature branch of the 
free energy per site / as obtained from thermodynamic integration, which yields J/kBT tr (D = 12) = 0.2511(10). The relative 
deviation |/wls — / ± |//wls between the thermodynamic integration and the WLS result is plotted in the inset. 



and 



B. Normalization of the DOS 



{0 if v ^ surface 

if v G surface, (5) 
/i(^)fzsurfacc 

This results in a number of M = 27 different classes. 
Within the context of N-fold way WLS the probability of 
any spin of a class i being overturned is then given by 



, A , n(C,AEi) 



,...,M, (6) 



where n(C, AEi) denotes the number of spins of config- 
uration C which belong to class i and pc—>c is given by 



Pc^c = 



min ( 1 9(Ec] 



if E c - G L 



sub 



if Ec <? i S ub, 



(7) 



where / su b denotes the considered energy subinterval over 
which the DOS is sampled and Ec — Eq + AEi. Classes 
are now chosen as follows. Firstly, one computes the 
integrated probabilities for a spin flip within the first m 
classes: 

Q m = ^P(A£;), m = l,...,M and Q = 0. (8) 

By generating a random number < r < Qm one 
then finds the class m from which to flip a spin via the 
condition Q m -i < r < Q rn . The spin to be overturned is 
chosen from this class with equal probabilities, whereby 
\og w g(E) and the energy histogram are now updated 
by means of the average lifetime r = 1/Qm- A detailed 
description of the algorithm was given in Ref. |16| . 



In order to estimate the DOS using WLS, the consid- 
ered energy range 



El JN G I = [E SIonad /JN, 0.2] , 



(9) 



where Aground = — [(3-D — 5) J + 4J S ]/JD, is the twofold 
degenerated ground state energy, was partitioned into 
several adjacent subintervals each containing an order of 
10 2 to 10 3 distinct energy levels, which were sampled on a 
Cray T3E in a parallel fashion using mostly 64 processors 
at a time. The DOS obtained from these simulations 
was then matched at the edges and suitably normalized, 
which we will describe in detail below. For the system 
thicknesses D = 8 and D — 12, as well as for the largest 
choices of L in case of D = 6 (L = 96, 128) only one run 
was performed over the entire energy range © denoted 
as basis-run, whereas all further runs have been restricted 
to a smaller energy range 



Ej JN G I c 



[E 1 /JN,E 2 /JN] 



(10) 



covering the mixed phase region in between the peaks 
of the doubly peaked energy distribution. As an illus- 
tration, / C ontor is marked in Fig. [S] by small arrows on 
the energy axis. Thus, the entire energy range @ is 
decomposed as I = 2i e ft U /center U / r i g ht , where we have 
/eft = [E gIound /JN,Ei/JN] and /right = [E 2 / JN, 0.2]. 
Correspondingly, one obtains the density of states g(E) 
by joining g(E) estimated for the intervals / c ft (taken 
from the basis run), /center, and / r i g ht (again taken from 
basis run). 

The single-spin-flip algorithm is more efficient in the re- 
gions covered by / cen ter which is due to the added ex- 
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L 


# runs 


In Zp = 






N 


N 


tazggf I 


6 


16 


6 


— 0.6932(5) 


-0.693147 


0.0076% 


6 


24 


3 


— 0.6931 (4) 


-0.693147 


0.0061% 


6 


32 


3 


-0.6932(3) 


-0.693147 


0.0082% 


6 
6 


48 
64 


3 
5 


-0.69318(5) 
-0.69311(6) 


-0.693147 
-0.693147 


0.0048% 
0.0049% 


6 
6 


96 
128 


2 
6 


-0.693112(2) 
-0.693144(9) 


-0.693147 
-0.693147 


0.0051% 
0.00042% 


8 


32 


4 


-0.6930(2) 


-0.693147 


0.027% 


8 


48 


2 


-0.69310(4) 


-0.693147 


0.0075% 


8 


64 


3 


-0.69312(6) 


-0.693147 


0.0038% 


8 


96 


1 


-0.69301 


-0.693147 


0.020% 


8 


128 


1 


-0.69304 


-0.693147 


0.016% 


12 


48 


3 


-0.69240(2) 


-0.693147 


0.107% 


12 


64 


1 


-0.692544 


-0.693147 


0.087% 


12 


96 


2 


-0.692686(3) 


-0.693147 


0.067% 


12 


128 


10 


-0.69281(4) 


-0.693147 


0.048% 



TABLE I: Logarithm of the partition function —(1/N) lnZp^o 
of a thin Ising film for different linear dimensions L and D, 
in case the density of states is normalized with respect to the 
ground state. The value in brackets states the standard de- 
viation. The exact value and the deviations from the latter 
are listed in the last two columns, respectively. For L — 32 
and D = 8 the run showing the largest deviation from the ex- 
act value (— (1/iV) In Z,3 = o = —0.692624) was excluded from 
data analysis. Then one has -(l/N)]nZ/3 =0 = -0.69307(9). 
Under # runs we have listed the number of independent sim- 
ulations. 



pense of the N-fold way algorithm concerning the book- 
keeping of classes. This was affirmed by a rough com- 
parison between both implementations for L = 128 and 
D = 12. The flatness parameter e varied between 0.8- 
0.95, and the final modification increment was usually 
of order Asg na i ~ 10~ 9 which yielded an overall simula- 
tional effort of order 10 6 MCS to 10 7 MCS for estimating 
the DOS over the range I©- As clear from the algorithm, 
WLS only yields a relative density of states, hence avail- 
able reference values must be employed in order to get 
the absolute DOS g(E). Normalizing the simulational 
outcome firstly with respect to the twofold degeneracy 
of the ground state, i.e., the free energy / will be exact 
for P — > oo, it instructive to examine how this accuracy 
for low temperatures carries over to infinite temperature 
(f3 — ► 0), where the partition function Z is dominated 
by the density of states around E = 0, and one has 
lima^o PF(0)/N = -{l/N)lnZ(0 = 0) = -In 2. Ta- 
ble 0] shows the latter quantity for all considered system 
sizes. As can be seen from the table there is an increasing 
deviation from the exact value with increasing width of 
the film D. While the results for D = 6 and D = 8 (the 
latter for small sizes L) agree with the expected value, a 
deviation for the larger system sizes, especially L = 128 
and D — 12 becomes apparent. This is related to a slow- 
ing down in the equilibration process in the multicanoni- 
cal (Wang-Landau) ensemble for decreasing modification 
increment Asi, as illustrated in Fig. which shows the 



visited states (E/JN, M/N) and the energy histogram 
H(E) recorded during Wang-Landau sampling for dif- 
ferent stages i of the simulation, where the modification 
increment Asi is used to modify the density of states. In 
case one has a small number of tunneling events during 
a certain simulation stage i, H(E) exhibits a kink at the 
barrier, since the stage is completed once the flatness cri- 
terion is fulfilled. Correspondingly, g{E) will suffer from 
large errors in the ordered phase, in case it is normalized 
using a reference in the disordered phase, and vice versa, 
errors will be enhanced in the disordered phase when us- 
ing a reference in the ordered phase (ground state). 
For D = 8 (excluding L — 32) and D = 12 we therefore 
employed the following approach. Utilizing the fact that 
one has at least random walk behavior for small modi- 
fication factor in each of the phases alone, we normalize 
the branch of the density of states g- (E) contributing to 
the low-energy ordered phase and the branch g + (E) con- 
tributing to the high-energy disordered phase separately, 
i.e., one has 



, g_(E) for E<E cut , 
<![L] 1 g + {E) for E /:,„ : . 



(11) 



where one obtains g~(E) by normalizing the simula- 
tional outcome g(E) with respect to the ground state 
g(E gmU nd) = 2 and is g+(E) obtained by normalizing 
g(E) with respect to the total number of states 



E-9(^) = 2 



L 2 D 



(12) 



In Eq. Ijllfl , E cut is taken to be the energy for which the 
energy probability distribution, estimated directly from 
the simulational outcome g(E), takes its minimum in be- 
tween the peaks at equal weight. Note, that in the sum 
Y,e9{E) = Y, E <E mi 9-(E) +T,e>e cM 9+(E) the term 
T, E <E cut 9-(E) is negligible. 

The additional error which is introduced by this normal- 
ization procedure then depends on the contribution of 
the mixed phase configurations to the energy distribu- 
tion (and the choice of E cut ). However, since these mixed 
phase configurations are exponentially suppressed at the 
transition point, the error is expected to be of the same 
order, and correspondingly the error due to the choice of 
_E cut as well. Note, that already for L = 32 and D = 8 
the double Gaussian approximation to the energy proba- 
bility distribution, which neglects any mixed phase con- 
tribution, provides a reasonably good approximation to 
the measured distribution, apart from small deviations 
in in the tails (cf. Fig.lBJ). 



C. Shape transitions 

In Ref. [l8|, Neuhaus and Hager addressed the se- 
vere problem of slowing down in simulations of first-order 
transitions in the multicanonical ensemble 19]. This was 
exemplified by studying the two-dimensional Ising model 
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FIG. 4: Visited states {E/JN, M/N) (left hand side) and 
energy histogram H(E) (right hand side), recorded during 
WLS in the space of energy over the interval E/JN € 
[—2.4414, —1.2207] for different stages i of the simulation. 
The simulation used 3.844 • 10 7 MCS in total, with a flat- 
ness criterion for the energy histogram of e = 0.9 and a final 
modification increment of Asfl na i ~ 4.139 ■ 10 -7 . 



on L x L square lattices (periodic boundary conditions) 
below the critical bulk temperature on the whole magne- 
tization interval [— L 2 ,L 2 ] whereby the sampling of con- 
figurations with magnetization M — Si was biased 
with the inverse probability distribution of the magne- 
tization g~ 1 (M). Specifically, it was found in Ref. [lg 
that these simulations suffered from a slowing down due 
to a discontinuous droplet-to-strip transition |2fj[. i.e., 
r oc exp(2i?crL), where t is the tunneling-time between 
droplet and strip configurations, a is the interfacial ten- 
sion, and R was measured to be R = 0.121(14). Note, 
that one has R ps 1 for non-multicanonical simulations. 
Of course, one needs a fairly good approximation to 
g(M), in order to sample the considered Hamiltonian in 
the multicanonical ensemble. Within the framework of 
WLS one may therefore simulate the system at a cer- 
tain inverse temperature /3 of interest by employing the 
flipping probability (single-spin-flip Metropolis) 



Pc^C 



1. 



g(M c ) 
g(M e >) 



exp(-/3 [Ea-Ee]) 



(13) 



for the transition from the state C to the state C Each 
time a state with magnetization M is visited, one updates 
g(M) according to g(M) — > g(M)- fi in complete analogy 
to the case where g(E) is used. Once this procedure has 
rendered g{M) accurate enough, one makes a production 
run, where g(M) is not altered anymore. Thermody- 
namic quantities can then be obtained by reweighting to 
the canonical ensemble. 

For the first order interface transitions in thin Ising-films 
as studied here, we have found evidence, that geometri- 




D=12 

FIG. 5: (b) Selected part of a time-series of the total magneti- 
zation per spin m — M/N as produced by WLS (single-spin- 
flip) in the space of magnetization, (a) Droplet at surface 
n — D where the positive field Hd/J = 0.25 acts, (c) Per- 
colated strip- like droplet. Note that in (a) and (c) only the 
positive spins are displayed as small spheres. Those spins 
closest to the shown L x L-surface are the lightest. 



cal transitions in the ensemble realized by Wang-Landau 
sampling in the space of magnetization, indeed hamper 
the simulations. While this poses no problem for the 
smaller systems like D = 8 and L = 32 where WLS us- 
ing g{M) yields very good results (cf. Fig.|5J), we observe 
pronounced effects for the largest considered system size. 
This is shown in Fig. Elb) where part of a time series 
is depicted which was recorded for D = 12 and L = 128 
during WLS sampling in the space of total magnetization. 
(Note also that the slowing down is more severe in case 
WLS using g(E) is employed for the same system size 
as obvious from Fig. 0]) The simulation was restricted 
to the interval m = M/N e [-0.55949,-0.45776] af- 
ter monitoring the time series of m for a much larger 
interval m 6 [—0.91553,0.10173] where Asj decreased 
from 5.0 • 10~ 3 to 7.629 • 10~ 8 over a simulation time 
of 1.632 ■ 10 7 MCS at J/k B T = 0.249719. The dis- 
tribution g(M) was then further iterated on the inter- 
val m G [—0.55949, —0.45776] where Asi was refined 
from 1.0 • 10~ 5 to 1.953 • 10~ 8 within 7.27 ■ 10 6 MCS 
and finally held fixed such that the depicted time se- 
ries could be recorded. Configurations were thereby 
monitored along the estimated position of the barrier 
m ps -103000/iV = -0.523885. Fig. GJa) and (c) show 
snapshots of the two possible coexisting structures which 
are the three-dimensional analogs (in the presence of a 
surface) to the droplet and strip shapes as studied in 
Ref. ■ In case Wang-Landau sampling is performed 
in energy space, the governing mechanism of the slowing 
down is not determined by now. From the joint energy- 
order parameter distribution (Fig. however, recorded 
for Wang-Landau sampling in the space of energy, one 
can at least conclude that one suffers from the fact that 
the ordered and the disordered phase are not distinctly 
separated in energy, as can be seen by inspecting the 
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P[n,mJ 




FIG. 6: Joint energy-order parameter distribution as obtained 
from WLS in the space of energy for a system size of D = 6 
and L = 16. The distribution was recorded using a fixed DOS 
g(E), which was taken from an usual adaptive WLS. 
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FIG. 7: Logarithm of the energy density of states logio(g(E)) 
for thicknesses D = 6, 8, 12 and linear dimensions L = 
48, 128. Smaller choices for L (in case of D — 6 and D = 8) 
are omitted in order to preserve clarity. Also indicated is the 
region where E cut , appearing in Eq. 1111 . is typically located. 
Here, both branches of the density of states gr_ and g+, are 
joined (D = 8, 12). In case of D = 6, g(E) was normalized 
solely with respect to the ground state degeneracy. 



spin 



") = 1^3) E EU 9( E ) e*P{-0E), (14) 



and the specific heat 



N 



«e 2 ) - (e) 



(15) 



Furthermore, important quantities like the free energy 
per spin can be directly computed 



-f3E 



(16) 



and the entropy per spin can be obtained from the inter- 
nal energy (|14|l and the free energy (|16fl 



f 



T 



(17) 



By measuring microcanonical averages (-)e during the 
last stage of a one-dimensional random walk in energy 
space, where g(E) is updated with the smallest increment 
Asfi na i we can also compute canonical averages of the 
order parameter (and higher moments) 



N 



(18) 



i.e. 



UH n )E9(E)e 



-0E 



,n\ E 



E 



(19) 



Thus quantities like the finite lattice susceptibility x 
N 



X 



((m 2 )-(|m|) 2 ), 



(20) 



as well as the fourth order cumulant U4 on which we 
concentrate in the following and which is defined as 



distribution of magnetization (along lines of constant en- 
ergy) which shows a noticeable three-peak structure for 
a range of energies e/J. Further studies are clearly nec- 
essary in order to clarify whether there are connections 
to droplet related phenomena. 



III. 



SIMULATION RESULTS 



A. Thermodynamic Quantities 

From the simulated DOS, as depicted in Fig. [7] we have 
calculated the first and second moment of the energy per 



U 4 = 1 



3{m 



2\2 : 



(21) 



become accessible. 

The distinctive feature of first-order phase transitions 
are phase coexistence and metastability. For the in- 
terface localization-delocalization transition considered 
here, this is reflected by jump discontinuities in the in- 
ternal energy (e) as well as the (absolute) magnetization 
(|m|) per site, which are depicted in Figs. I§1 and UHl 
respectively, and also by hysteresis effects encountered 
when heating and cooling the system as exemplified in 
Fig. |3J Considering the internal energy (Fig. |5J| for fixed 
D and varying linear dimension L, one can clearly see 
that one actually does not observe discontinuous jumps 
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e/J 



FIG. 8: Panel (a) shows the double Gaussian approximation l|25fl to the energy probability distribution Pl,d(c) for the system 
of size L — 32 and D — 8 at the finite volume transition point ((3tr(L,D) = 0.247255(10)) as obtained from WLS in the space 
of magnetization (single-spin-flip) by reweighting back to the canonical ensemble. Panel (b) shows the corresponding full joint 
energy- order parameter distribution Pi,D(e, rre), while (c) shows the projection onto the magnetization axis. 
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FIG. 9: Internal energy (e) for different linear dimensions L and film thicknesses D. Estimates for the inverse temperature 
J/3tr(D) — J/kBT tT (D) of the triple point are indicated by arrows. The horizontal solid lines mark the value (e + + 2e~)/3 
where the curves are expected to cross. In (b) the data obtained from WLS using g(M) (D = 8 and L = 32) is plotted for 
comparison. Here, g(E) for D = 8 and L = 32 was normalized solely with respect to the ground state. Within the inverse 
temperature range displayed in the inset of (a), the average relative errors in (e) for D = 6 amount to 0.17%, 0.54%, 0.55%, 
0.42%, and 0.45% for L — 16, 128, respectively. For D = 8 and L = 32 the average error amounts to 0.18% in the range 
0.2455 < J/k B T < 0.2485, when using g(M) while it is 1.0% within the same range when using g(E). Note, that for D = 8 
and L = 96, 128, as well as for D = 12 and L — 64 the DOS was determined only once, hence no error bars are displayed. 



of the quantities in question but a continuous behavior 
that sharpens to the asserted step-like behavior with in- 
creasing linear system size L. This rounding is related to 
the fact that a true phase transition can only occur in the 
thermodynamic limit, where in equilibrium, approaching 
the transition temperature from above, the energy of the 
system discontinuously jumps from e + (interface in the 
center of the film) to e~ (interface tightly bound to the 
wall), while for a finite volume the system may jump back 
and forth between the latter states and the observed equi- 
librium behavior is thus continuous in temperature. The 
rounding of the transition in finite systems can also be 
observed for the specific heat c depicted in Fig. ITT1 which 



exhibits narrow peaks that are remnants of the ^-function 
singularities one would get when differentiating the dis- 
continuous energy in the infinite volume limit. Apart 
from the finite size rounding, one can see that the posi- 
tions of the maxima of c and the minimum of the fourth 
order cumulant U4 (Fig. I12fl are systematically shifted 
towards higher /?- values for increasing linear dimension 
L. 

From the crossings of the energy curves for different linear 
dimensions L, one can get a first idea about the achieved 
accuracy for the different film thicknesses D, because 
they should cross to a very good approximation in the 
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FIG. 10: Average absolute magnetization per spin (|m|) of 
a thin Ising film for different linear dimensions L and film 
thicknesses D. Within the inverse temperature range dis- 
played in the inset of (a), the average relative errors in (|m|) 
for'D = 6 amount to 1.0%, 5.2%, 6.3%, 6.1%, and 9.3% for 
L — 16, 128, respectively. 



point 21 



(p tT (D), (e+ + 2e-)/3), 



(22) 



where Pti(D) is the infinite system transition point. 
Hence, the crossing points for different L 



(e(f3 CIOSa ,L,D)) = (e(/3 

cross • 



(23) 



actually provide an estimator for the infinite system tran- 
sition temperature, which is expected to deviate from 
f3tr(D) only by an amount exponentially small in system 
size |21| ■ As can be seen from the inset of Fig. Ha) in 
case of D = 6, the various crossings are indeed scattered 
in a narrow region around the extrapolated infinite vol- 
ume transition point for L > 32. For smaller values of L 
exponential corrections still make a noticeable contribu- 
tion. For the larger thicknesses D > 8 the region where 
the energy curves cross is noticeably larger. Particularly, 
one obtains that the errors resulting from averaging over 
different runs are too small to fully account for the devi- 
ations (excluding L — 32 for which g(M) was employed). 
This is related to the fact that for the thicknesses D = 8 
and D = 12 only a single run was performed over the 
entire energy range 10 while further runs were restricted 
to the mixed phase region in between the peaks, because 
the slowing down, as described in the preceding subsec- 
tion, was not foreseen. When one uses the normalization 
condition jlip. the proper strategy would certainly be to 
enhance the simulational effort in the pure phases, down 
to the ground state and up to E = 0, since the reference 
density of states is known for T = and f3 = 0. This 
is necessary, in order to minimize the accumulation of 
errors in the density of states, since the Wang-Landau 



method and similar adaptive algorithms, do in general 
not exhibit an error distribution that is flat in energy |36| . 
Hence, for D — 8 (excluding the simulation using g(M)) 
and D = 12, we believe the true errors to be larger than 
the error bars displayed in Figs.l9Tb')-fc).[T2Tb). [TlT b1-fc1 
and when quantitatively referring to errors of the ther- 
modynamic quantities, we thus restrict ourselves here to 
D = 6, where we have reliable error estimates. 
Exponential corrections to the crossing points are pre- 
sumably much smaller than the scatter in the energy 
crossings for D > 8 and one may therefore conclude that 
the deviations in the crossings for D > 8 are not due to 
corrections to scaling, but reveal the actual error in the 
density of states for this region. This is also the case 
for the other quantities like the specific heat for example 
(Fig. Ilir bVfc')'). Thus, the analysis of the systems with 
larger thicknesses D = 8 and D = 12 is certainly more 
difficult and less accurate. 

One can however roughly estimate the order of magni- 
tude of the latter uncontrolled error, which also serves to 
support the above picture. For example, from the density 
of states of the largest system (D — 12 and L = 128), 
one can estimate, that a relative error in the density of 
states g(E) of the order ~ 10 _1 (referring to the results 
for the 50 x 50 2D Ising model in Ref. [16| this seems to 
be a reasonable assumption) , in the narrow region corre- 
sponding to the peak of the ordered phase of the energy 
probability distribution, can result in a displacement A/? 
of the peak position /3 c max of the specific heat and also of 
the step location of the internal energy, which is approxi- 
mately of the order A/3//3 c max ~ 10 -4 . In case of D = 12 
and L = 48, a relative deviation of this order could al- 
ready be caused by a relative error in g(E) which is of the 
order ~ 10~ 2 in the above region. These considerations 
comply well with the observed scatter. 



B. Finite size scaling 

When one deals with second order phase transitions, 
the characteristic feature is a divergent spatial correla- 
tion length £ at the transition point (3 C (where one ob- 
serves fluctuations on all length scales) implying power- 
law singularities in thermodynamic functions such as the 
correlation length, magnetization, specific heat and sus- 
ceptibility. This is in sharp contrast to a first order tran- 
sition where the correlation length in the coexisting pure 
phases remains finite and concerning finite size scaling 
the volume of the system turns out to be the relevant 
quantity. For a thin film geometry where one has fixed 
D and varying linear dimension L, finite size scaling will 
thus involve the quantity L 2 . This can be shown by ap- 
proximating the energy distribution Pl,d{z) of the pure 
phases by a Gaussian 0, 122I |23| centered around the 
infinite-lattice energy per spin (e) 



Pl,d (g) 



' L 2 D 
2Trk B T 2 c 



exp 



2k B T 2 c 



L 2 D 



(24) 
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FIG. 11: Specific heat c of a thin Ising film for different linear dimensions L and film thicknesses D. In the interval 0.2400 < 
J/k B T < 0.2415 the average relative error for D = 6 amounts to 3.9%, 9.6%, 12.0%, 5.4%, and 15.9% for L = 32, 128, 
respectively. For D = 8 and L = 32, Wang-Landau sampling in g(M) yields an average error of 2.6% within the range 
0.2455 < J/k B T < 0.2480. Note, that we have no statistics for D = 8 and L = 96, 128 (b), as well as for D = 12 in case of 
L = 64 (c). 




FIG. 12: Reduced fourth order cumulant of a thin Ising film for different linear dimensions L and film thicknesses D. The inset 
of panel (a) shows the region where the cumulants for the various linear sizes L cross (D = 6). In the vicinity of the minima 
positions, the relative errors in Ui in case of D = 6 amount to 11%, 30%, 17%, 11%, and 7% for L = 24, 128, respectively. 
WLS in the space of total magnetization M with fixed g(M) yields an error of 4% (D = 8, L = 12). Note, that for the data 
corresponding to D = 8, as plotted in (b), we have no statistics for L = 96 and L — 128, i.e., for the latter sizes the DOS was 
estimated only once. 



where c denotes the infinite-lattice specific heat. Since 
one has phase coexistence at a first-order transition, 
the probability distribution of the energy will be dou- 
ble peaked at the transition point /3t T (D) = l/kBT tI (D), 
where (e) jumps from e~ (low energy phases, interface 
at one of the two walls) to e + (single high energy phase, 
interface centered in the middle of the film), i.e., the free 
energy branches / intersect at a finite angle in the in- 
finite system, as can be seen from Fig. I13f a). when in- 
specting the curves around the transition point (cf. also 



Fig 3(b)) It is essentially this non-analyticity in the free 
energy, which gives rise to the discontinuous behavior 
of the internal energy. In a finite system however, the 
free energy remains differentiable and the intersection is 
rounded. 

Hence, at the transition point, Pl,d{g) is a superposi- 
tion of two Gaussians l|24|l centered at (e) = e , while 
slightly away from the transition at T — T tr + AT they 
are centered at energies + c^AT, where are the 
specific heats in the disordered (+) and ordered phases 
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(-), which are assumed to be constant in the vicinity of 
the transition, i.e., for sufficiently small AT. Each of the 



J 



Gaussians is then weighted by Boltzmann factors of the 
corresponding free energies f^, and one thus arrives at 



Pl,d{z) = A 



exp 



[e - (e+ + c+AT)] 5 
2fc B T 2 c+ 



L 2 D 



r 



exp 



-AT)] 5 



2k B T 2 c- 



L D 



(25) 



where the weights a are given by 

?+ 



q exp 



=p- — L 2 -D 

2fc B T 



and ^4 reads 



.4 



exp 



(/ + + n 

2fc B T 



Z/T> 



/ L 2 D 
2irk B T 2 ' 



(26) 



(27) 



Since we have a single high energy phase and two low 
energy ordered phases we set q + — 1 and q~ = q = 2 
in the following. At the transition all phases have equal 
weight [2lll2 1| such that the area under the peak at e~ is 
q times the area under the peak at e + which is satisfied 
by Eq. Within the framework of the Ansatz 

one then proceeds by calculating the energy moments as 
usual via 22] 



Jde'e' n P L , D (e') 



(28) 



Computing then (e) 
Eq. J5S|) we obtain 



Jde'P L , D (e') 
at the transition point by means of 



qe 



1 



(29) 



(horizontal lines in Fig. which is exact, apart from 
exponential corrections due to mixed phase contribu- 
tions which are neglected in the double Gaussian ap- 
proximation. Upon using the fluctuation relation 1151) 
orc = d(e) /dT in conjunction with Eq. (|28() one can cal- 
culate the specific heat to leading order 



2 + + a c 



- a 
e~ 



+ (c+-c-)AT] 2 aV 



(a+ + a-) 2 



kp,T 2 



L 2 D, (30) 



which is seen to take its maximum for a + = a~ in 
Eq. I|25(l . The position of the latter is thereby shifted 
away from the infinite lattices transition temperature by 
an amount of 



AT = T C ^(D, L) ~ T tr (D) = k B T 2 . 



2 m <? 



1 



AeD L 2 ' 

and the height of the peak is found to be 
c+ + c~ Ae 2 T> ~ 



4fc B T 2 



(31) 



(32) 



where Ae = e + — e~ is the latent heat. For convenience 
we may re-express Eq. (),'.'> 1^ in terms of the inverse tem- 
perature (3 = l/fc B T which yields 



c ™*(D,L)=(3 tI (D) 



ln2 1 
AeDL 2 ' 



(33) 



Thus, the inverse temperature /3 c max (D) at which the spe- 
cific heat peaks, provides a definition for a finite-lattice 
(pseudo) transition temperature from which the infinite- 
lattice transition temperature can be estimated via finite 
size scaling, i.e., by extrapolating L — > 00. 

A similar argumentation ap plie s to the distribution of 
the order parameter Pl i _d(to) [23 yielding the same scal- 
ing behavior for the susceptibility \i i- e - 



p x ^(D,L)-p tr (D) cx (tfD)- 1 



X' 



cx L 2 D, 



(34) 
(35) 



and one can show that the fourth order cumulant U4 
(Fig.UU takes a minimum value at an inverse tempera- 
ture /3;ymm(Z?, L) which is again shifted 



/3™.„(A£)-Ar(£)cx (L 2 D)~ 



while the minimum [/I 11111 



obeys 



-L 2 D. 



(36) 



(37) 



Furthermore it was shown [2£j , that the shift in the cross- 
ing points of the cumulants for different system sizes is 
proportional to N~ 2 , which is negligibly small on the 
scale of N = L 2 D. 

Fig. I14f b) now shows the maximum values of the re- 
sponse functions c max , x maX ) an d the minimum U™ m of 
the cumulant as function of 1/L 2 for the three differ- 
ent thicknesses D = 6, 8, 12. As can be seen from the 
plots, the data comply well with the behavior predicted 
by expressions (|32f) , (|35|l and l|37|) . Considering the fourth 
order cumulant U4 in case of I? = 6, one observes that 
sub-leading corrections to scaling are still present for the 
smaller linear dimensions T, but the expected linear be- 
havior in L 2 is born out for the three largest choices of 
L. 

The definition for the finite lattice transition tempera- 
ture considered so far, e.g. Eq. involve leading order 
corrections of 1/L 2 . An alternative definition of the tran- 
sition temperature which has the additional benefit that 
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FIG. 13: (a) Free energy per spin / of a thin Ising film for different linear dimensions L and film thicknesses D. Note, that in 
(a) only the data for L = 128 is plotted, while (b) shows / on a finer scale. In the range 0.23 < J/ksT < 0.27, the average 
relative error in the free energy (D = 6) is 0.0143%, 0.011%, and 0.00026% for L = 16, 32, and 128, respectively, (c) entropy 
per spin s for D — 6. The error in s amounts to 0.64%, 0.64%, 0.47%, and 0.47%, within the range depicted in the inset of 
panel (c). 




FIG. 14: Panel (a): extrapolation of peak positions /3 m ax,min(-D, L) of the specific heat c max and the fourth order cumulant 
Uf 1111 for the different film thicknesses D. Panel (b): maxima of specific heat c max and susceptibility x maX i as weu as minimum 
of the fourth order cumulant U4 1 ' 11 as function of L 2 . 



the latter corrections are absent was given in Ref. [26(. 
Here, it is utilized that at the infinite-lattice transition 
point Ptr(D) all phases coexist which implies that the 
sum of the weights of the q ordered phases equals q times 
the weight of the disordered phase, i.e., 



K(p ev! ,L,D) = — — - — -— q, (68) 

Le>e cut PL,D(e,(i cw ) 



where -Pl, 23(e) is the (finite-size) energy probability dis- 
tribution, and f3 evr (D, L) differs from Ptr{D) only by cor- 
rections exponentially small in system size. The energy 
e cut appearing in Eq. I|38|) is taken to be the internal en- 
ergy at the temperature where the specific heat is maxi- 
mal m. 



C. Transition temperatures 

Now, we can extract the infinite volume transition 
point Pti(D) from the finite size data, i.e., as Eqs. 
and (|36fl suggests by fitting the peak positions for fixed 
D to 

/3max,mm(A L) = /3max,min(A Oo) + j^, (39) 

where /3 m ax,min(-Dj L) stands for the location of the max- 
imum of the specific heat /3 c ma X (D,L) and the location 
of the minimum /3jjmm(_D, L) of the fourth order cumu- 
lant at finite L, while /3 m ax,min(-D, 00) denotes the infi- 
nite volume limit (L — > 00) of the corresponding inverse 
temperatures, which is an estimate of the infinite sys- 
tem transition point @ tT (D). Alternatively, we have also 
employed the finite volume estimator /3 ew (D,L) of the 
transition point, as defined by the condition l|38|) . 
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D Ae{D)/J J/3 c max(D) J/3 J/3 ew (D) JA^D) 



6 0.261(6) 
8 0.300(2) 
12 0.236(3) 



0.24082(2) 0.24079(4) 
0.24726(3) 0.24716(7) 
0.25115(4) 0.25109(4) 



0.24082(1) 0.24082(4) 
0.24725(2) 0.24725(10) 
0.25117(2) 0.25115(10) 



TABLE II: Estimates for the latent heats Ae(D) and the 
inverse transition temperatures of the first order interface 
localization-delocalization transition for different film thick- 
nesses D. /3 c max(_D, oo) and f3 um i n (D, oo) are the estimates of 
the transition point /3 tr (D) originating from an extrapolation 
of peak positions as described in the text, while /3 ew (D, oo) 
denotes the estimate from the equal weig ht rule (EHJ. The 
final estimate of the inverse temperature {3 t r(D) of the triple 
point is stated in the last column. 



The individual results for the infinite system transition 
points are summarized in Table ITT1 In the last column of 
Table ITT1 we state our final estimate of the infinite system 
transition point /3 t r(D), based on weighted averages over 
the estimates listed in column 2-4. Concerning the error 
in our final estimate of fitr{D) we have also accounted 
for the scatter in the crossings of the energy curves as 
depicted in Fig. [5| and the crossings in the fourth order 
cumulant U4, see Fig. ^] While we find that the or- 
der of magnitude of the error as determined from the 
various finite lattice estimators considered above, com- 
plies well with all the data for D = 6, especially the 
latter crossing points, we may have uncontrolled errors 
in case of the larger thicknesses D = 8 and D = 12, due 
to the aforementioned lack of statistics deep in the pure 
phases. In these cases we consider here as a conservative 
error estimate the extremal crossing points as an upper 
bound to the transition point, which results in the error 
of ftu (D > 8) as given in the last column of Table ITT1 
Fitting the locations of the maxima of the specific heat 
to Eq. as depicted in Fig. lTlT a). one can also deter- 
mine the latent heat Ae which is however less accurate 
than computing Ae from the distribution Pj, d(E) via 

M 



Ae(L, D) = Ae(.D) + const x L~ 



(40) 



which yields the values stated in column 2 of Table ITT1 
Concerning the extrapolation l|39|) of the positions of the 
minima {7™ m and the maxima c max we have used only 
data for L > 32 in case of D = 6. For these lattices, 
exponential corrections to (3 ew (D,L) cannot be resolved 
within the achieved accuracy. This is also the case for the 
larger film thicknesses D and all choices of L. Hence, the 
values listed in Table HT1 for [3 CW (D) are simply averages 
over the various lateral system sizes L (L > 32 in case of 
D = 6). 



D. 



In 



Wetting Temperature of the semi-infinite 
system 



order to determine the wetting temperature 
/9 W = liniD^oo (3tr(D) of the semi- infinite system, we have 
studied Hamiltonian with D = 12 and L = 48 along 
the branch of positive bulk-magnetization at the inverse 
temperature f3 — 0.251 near the expected location of the 
wetting temperature /3 w (JZi). We have performed sim- 
ulations for five different sets of surface fields (symmet- 
ric, i.e., Hi = H D ), namely H x /J = -0.25, -0.125, 0, 
0.125, and 0.25, utilizing a conventional Metropolis al- 
gorithm in order to measure the surface magnetization 
= (E ies urfacci St)/N using up to 10 7 MCS for av- 
eraging. This selection of surface fields allows one to 
reweight to all fields in the range [—0.25 J, 0.25 J] for a 
range of inverse temperatures J/3 <E [0.249,0.253]. Note, 
that the metastability is strong enough (cf. Fig. |3J) that 
the system remains in the ordered phase (initially all 
spins up) even for Hi/ J = —0.25. According to the 
Young equation 29] the walls are wet by spin down, if 
the difference Aer w between the surface free energy of the 
wall with respect to a positively magnetized bulk a w + 
and the surface free energy against a negatively magne- 
tized bulk cr w _ exceeds the interfacial tension a of the 
3D Ising-model |2Sj at an infinite distance from the wall. 



Act. 



•>w+ 



> CT 



(41) 



By symmetry ct w _(— Hi) equals a w+ (Hi), i.e., the free 
energy cost of a wall favoring spin up with respect to 
a positively magnetized bulk. Thus we can perform a 
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FIG. 15: (a) Shown are the interfacial tension a of the 3D 
Ising-model (HP) taken from Ref. @, fitted by an 8 th degree 
polynomial in order to smoothly interpolate between the data 
points as well as the the quantity Acr w appearing in the Young 
equation H4H . The position of the crossing point yields the 
wetting temperature J/3 w (-ffi) = J/fc B Tw(iii) = 0.25212(5). 



14 



thermodynamic integration 30] 

A<t w = a w+ (-Hi) - cr w+ (ifi) 



ffi 



dff^miGffO)^ i?i=0.25J (42) 



and determine the wetting temperature (3 w (Hi) by the 
condition Aa w — a, which yields J(3 w (Hi) — 0.25212(5) 
as depicted in Fig. HTT a). 

Describing the semi-infinite system by means of the wet- 
ting film thickness / leads to the effective interface po- 
tential 

V e g(l) = aexp(— k/) — 6exp(— 2nl) + cexp(— 3k/), (43) 

which has the meaning of a free energy cost when plac- 
ing a (flat) interface at distance I from the wall. Upon 
minimizing V e g(l) with respect to I one finds the equi- 
librium position of the interface. Eq. (|43H includes only 
the lowest powers of exp(— nl) which are necessary to de- 
scribe a first order wetting transition in the semi-infinite 
system. The coefficient a explicitly depends on temper- 
ature, while the temperature dependence of b and c is 
neglected (c > in the following) [37| . All coefficients 
have the same magnitude as the interfacial tension be- 
tween bulk phases and one finds a first order wetting 
transition for b > at a w = b 2 /4c, where the interface 
jumps discontinuously into the bulk 0,1^]. Now, for a 
film one has an additional contribution from the second 
wall and the effective potential reads [32| 



AKff,Fiim(0 = V cS (l) + V cS (D-l)-2V cS (D/2) 

™2/~2 a2 , . -21 



with 



and 



= c [to 2 (to z — r) z + tfhf' 

_b- 6cexp(-«D/2) 
~ 2~c ' 

a — a w — &exp(— nD/2) 



(44) 



(45) 



(46) 



In Eq. I|44(l we have utilized the auxiliary variable 



m = 2exp(-K£>/2){cosh[K(Z - D/2)] - 1} 
= (cxp(-KD/4)n[l- D/2}) 2 

+higher orders of [I - D/2] . (47) 

In the film, r > gives rise to first order interface 
localization-delocalization transitions and t = then de- 
notes the triple temperature. Hence, for large D we have 
from Eq. J32J) 



atr = a wc t + bexp(~nD/2), 



(48) 



i.e., the triple temperature differs from the wetting tem- 
perature only by a term exponentially small in kD/2 and 
is larger than the wetting temperature (b > 0). Within 



mean field theory k would have to be identified with the 
inverse bulk correlation length £b However, from the 
two-field Hamiltonian approach developed in Ref. [^11 we 
know that k/2 has to be replaced by 



2£b0' 



= 1 + uj cS /2, 



(49) 



where w c ff is the effective wetting parameter which be- 
comes lim T ^ T + w, 



: fF = k^T/A-Ka^ upon lowering the 



temperature T towards the wetting temperature T w |34j . 
From a simple exponential fit of the form l|48|l we get 
k/2 = 0.430(8). (Note, that this has to be regarded as 
an effective value since we neglect any temperature de- 
pendence of k within our range of triple temperatures 
fi tr (D)). Evaluatin g no w at T w /T cb — 0.88 where we 
employ £b ~ 0.88 [28|. yields 9 ~ 1.3, which is com- 
patible with the values extracted for 8 by Parry et al. 
|34| and clearly differs from the value 9=1 expected 
from mean-field theory. Of course, making more quanti- 
tative statements would require data from additional film 
thicknesses D, but the above considerations clearly indi- 
cate that our data nicely supports the asserted functional 
dependence of Pti(D) on D, i.e., Eq. 



IV. 



CONCLUSION 



We have studied the interface localization- 
delocalization transition in a thin Ising-film JQ) for 
a choice of parameters, where the transition is pro- 
nounced first order for all studied thicknesses D = 6, 
8, and 12. Checking for the correct behavior of the 
logarithm of the partition function In Z which should 
converge to iVln2 as (3 — > 0, we find reasonable 
agreement for D = 6 within error bars (cf. Tabled) 
In contrast, for D > 6 we see rather clear deviations 
from the expected value with relative deviations up to 
10~ 3 . We attribute this behavior to a slowing down 
encountered in the flat energy-histogram ensemble. Dif- 
ficulties also arise, when one considers to sample a flat 
magnetization distribution, although simulation results 
suggest that the slowing down is less severe. Here, we 
find evidence for a discontinuous shape transition, as 
studied by Neuhaus and Hager F° r t ne larger 

thicknesses (D > 6) we therefore suggest to employ 
an additional reference for the disordered phase (total 
number of states), in order to get the proper relative 
weight between the coexisting phases, thus correcting 
for the lack of tunneling events, in the late stages of 
the algorithm. The triple temperatures (3t r (D) of the 
interface localization-delocalization transition can then 
be determined with a relative accuracy of the order 
10~ 4 while the relative error in the latent heats is of 
the order 10~ 2 . The triple temperatures are seen to 
differ from the wetting temperature of the semi-infinite 
system by a term exponentially small in film thickness 
D as predicted by the sharp-kink approximation to the 
capillary wave Hamiltonian, provided the length scale 
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n is identified with the results of Parry and co-workers, 
i.e., Eq. g5J. 

When one compares the present results based on Wang- 
Landau sampling [l3L[T^.IT5llT^ | to the first study of first 
order interface localization-delocalization transitions Q 
where simple Metropolis and heatbath Monte Carlo 
algorithms were used, a major improvement of accuracy 
is clearly seen. On the other hand, the systematic prob- 
lems due to entropic barriers described in our work show 
that it would be problematic to apply the Wang-Landau 
algorithm to larger systems than used here. Note, that 
the largest sizes used by us, 128 x 128 x 12 - 1.97 • 10 5 
Ising spins, distinctly exceed the sizes analyzed in most 



previous applications of this algorithm [Lj, [LJ, Il5l [1 
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